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Parameterisation  of  eddies  in  coarse  resolution  models 

Peter  D.  Killworth 

Southampton  Oceanography  Centre,  Southampton,  England 


Abstract.  Some  requirements  on  eddy  parameterisations  are  discussed, 
especially  the  implications  of  expressing  them  in  terms  of  the  quasi-Stokes 
velocity  and  the  modified  mean,  rather  than  Eulerian  mean,  density.  The 
difference  between  the  two  means  is  second-order  in  perturbation  amplitude 
and  thus  small  in  the  fluid  interior  (where  formulae  to  connect  the  two  exist). 
Near  horizontal  boundaries,  the  differences  become  first  order,  and  so  more 
severe.  Existing  formulae  for  quasi-Stokes  velocities  and  streamfunction 
also  break  down  here.  The  layer  in  which  the  largest  differences  between 
the  two  mean  densities  occur  is  the  vertical  excursion  of  a mean  isopycnal 
across  a deformation  radius,  at  most  about  20  m thick.  Most  climate  models 
would  have  difficulty  in  resolving  such  a layer.  It  is  shown  that  extant 
parameterisations  appear  to  reproduce  the  Eulerian,  and  not  modified  mean, 
density  field  and  so  do  not  yield  a narrow  layer  at  surface  and  floor  either. 
Both  these  features  make  the  quasi-Stokes  streamfunction  appear  to  be 
non-zero  right  up  to  rigid  boundaries,  so  that  we  must  query  what  are  the 
relevant  surface  and  floor  quasi-Stokes  streamfunction  conditions,  and  what 
are  their  effects  on  the  density  fields.  To  answer  this,  a variety  of  eddy 
parameterisations  is  employed  for  a channel  problem,  and  the  time-mean 
density  is  compared  with  that  from  an  eddy-resolving  calculation.  The 
parameterisations  were  only  successful  if  the  vertical  component  of  the  quasi- 
Stokes  velocity  vanished  at  top  and  bottom  as  in  current  practice,  but  all 
were  almost  equally  successful  given  proper  tuning.  One  parameterisation, 
based  on  linear  instability  theory,  is  extended  to  a global  geometry.  In  low 
and  mid-latitudes,  because  the  predominant  orientation  of  the  instability 
wavevector  is  north-south,  the  main  quasi-Stokes  flow  is  east-west,  only 
becoming  the  more  traditional  north-south  at  higher  latitudes. 


1.  Introduction 

The  ocean  component  of  climate  models  is  necessar- 
ily of  coarse  resolution,  and  it  does  not  possess  eddies. 
This  paper  discusses  various  aspects  of  the  issue  of  rep- 
resenting eddies  in  such  models.  First,  it  discusses  the 
extant  parameterisations  (Section  2).  Section  3 exam- 
ines some  hitherto  unappreciated  details,  with  reference 
to  the  ‘modified  mean  density’  concept;  it  is  shown  that 
the  difference  between  this  density  and  the  Eulerian 
mean  is  largest  at  surface  and  floor.  The  region  that 
this  difference  occupies  is  too  shallow  to  be  resolved  by 
coarse  models,  so  that  eddy  effects  occur  in  what  ap- 
pears to  such  models  as  a delta-function.  Current  pa- 
rameterisations do  not  generate  solutions  which  include 
these  differences,  even  with  adequate  resolution.  Some 
eddy-resolving  channel  experiments  are  used  in  Section 
5 to  explore  whether  the  relevant  boundary  condition  is 
that  of  no  eddy  flux  (formally  correct  beyond  the  thin 
layer)  or  not;  the  conclusion  is  that  the  usual  boundary 


condition  must  be  employed.  Last,  a parameterisation 
based  on  linear  instability  theory  (introduced  in  Section 
4)  is  extended  in  Section  6 from  channel  simulations  to 
the  global  domain  and  is  briefly  compared  with  the  Gent 
and  McWilliams  (1990)  parameterisation.  The  direc- 
tionality of  linear  instability  at  low  and  mid-latitudes 
is  such  that  the  ‘bolus’  fluxes  are  oriented  east-west, 
rather  than  north-south,  at  these  latitudes,  because  the 
beta  effect  constrains  maximal  instability  to  a north- 
south  orientation  (which  yields  east- west  fluxes). 

2.  Background 

A variety  of  schemes  has  been  suggested  to  include 
eddy  effects  in  coarse-resolution  ocean  models.  These 
schemes  divide  into  two  categories.  The  first,  which  we 
shall  be  examining  here,  involves  adding  terms  to  rep- 
resent the  additional  thickness  flux  by  baroclinic  eddies 
( Gent  and  McWilliams,  1990;  Greatbatch  and  Lamb, 
1990;  Gent  et  al.,  1995;  Visbeck  et  al.,  1997;  Treguier 
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et  al.,  1997;  Killworth , 1997,  1998;  Greatbatch,  1998). 
The  second  (Neptune)  involves  a representation  of  the 
statistical  properties  of  eddies  on  the  mean  flow  ( Eby 
and  Holloway,  1994;  Merryfield  and  Holloway,  1997), 
and  will  not  be  discussed  in  detail  here.  Methods  to 
represent  propagating  features  (e.g.  Agulhas  rings)  do 
not  seem  yet  to  be  available. 

Eddy  parameterisations  have  been  designed  with  var- 
ious criteria  in  mind  by  different  authors.  Initial  re- 
quirements include  mimicking  both  baroclinic  instabil- 
ity - “reducing  APE”,  in  some  sense,  and  barotropic 
instability  - “reducing  lateral  shear” , as  well  as  the  Nep- 
tune approach  of  making  the  ocean  resemble  a rectified 
eddying  ocean.  All  these  are  quasi-steady  in  nature 
and  ignore  the  fact  that  instability  generates  variability 
on  weeks  to  decades  through  nonlinear  interactions.  I 
am  not  aware  of  any  parameterisation  which  attempts 
to  put  temporal  variability  into  the  ocean  (by  stochas- 
tic forcing,  perhaps),  so  that  coupled  ocean-atmosphere 
models,  even  with  parameterisations,  are  unlikely  to 
possess  a realistic  degree  of  self-induced  oceanic  vari- 
ability. 

Two  other  pragmatic  requirements  are  that  parame- 
terisations should  not  induce  erroneous  ocean  behaviour 
(e.g.  breaking  conservation  laws  such  as  momentum), 
and  should  ensure  that  the  numerics  of  the  models  be- 
have (the  original  concept  behind  eddy  diffusivity,  of 
course). 

The  manner  in  which  a parameterisation  is  couched 
depends  on  what  belief  structure  about  eddy  behaviour 
is  used.  In  the  literature  already  are  suggestions  for 
thickness  smoothing,  potential  vorticity  (q)  conserva- 
tion, energy  loss,  energy  conservation,  and  q smoothing. 
These  effects  are  usually  placed  in  the  tracer  equations, 
but  it  is  possible  to  include  effects  in  momentum  equa- 
tions also.  There  are  potential  difficulties  near  horizon- 
tal surfaces,  discussed  later. 

2.1  Formulations 

Most  authors  seek  an  “equivalent  to  the  bolus  veloc- 
ity” , namely  some 

(u+,u+,m+) 

which  is  fully  three-dimensional  and  non-divergent,  such 
that  u+u+  advects  (a  form  of)  the  density  adiabatically 
- i.e.  build  in  our  belief  that  density  (or  neutral  den- 
sity) effects  are  adiabatic.  This  is  usually  represented 
as  a 2-dimensional  streamfunction 

® = (V'l.lfe) 

such  that 

u+  = ipn,  v+  = if)2z , w+  = -Vff  • 'P. 


The  most  logical  approach  to  date  is  the  transient- 
residual-mean  (TRM)  theory  introduced  by  McDougall 
(1998,  and  earlier  references  therein);  McDougall  and 
McIntosh  (2001,  hereafter  MM)  give  more  detail  on  the 
same  material.  Another,  highly  related,  approach  is  to 
use  density-weighted  averaging  (cf.  Greatbatch,  submit- 
ted ms;  de  Szoeke  and  Bennett,  1993).  The  TRM  theory 
applies  to  low-pass  temporally  averaged  quantities  and 
deduces  a quasi-Stokes  velocity  u+  which  is  related,  but 
not  identical,  to  the  bolus  velocity.  (The  two  are  not 
identical  because  the  background  mean  flow  involves 
averages  on  two  different  surfaces,  though  they  are  fre- 
quently similar.)  Formulae  have  been  derived  for  small 
perturbations  by  McDougall  (1998)  and  MM,  involving 
only  averages  at  constant  depth.  The  quasi-Stokes  vec- 
tor streamfunction  is  given  to  second  order  in  amplitude 
by 

$ = + ^ (l) 

Pz  Pz  \PzJ  ’ 

where  the  suffix  II  denotes  the  horizontal  component, 
and  <j>  = (1/2 )p'2.  The  vertical  derivative  of  'P  is  the 
horizontal  component  of  u+ . The  second  term  is  usually 
small  compared  with  the  first  and  is  neglected  hence- 
forth. 

Since  eddying  motions  are  believed  to  conserve  den- 
sity, this  implies  that  the  definition  of  density  must  be 
modified.  McDougall  (1998)  shows  that  rather  than  us- 
ing the  Eulerian  mean  density  p at  a (vertical)  point 
(EMD  for  short),  one  should  interpret  density  as  be- 
ing the  inversion  of  the  mean  depth  of  a given  density 
(termed  the  ‘modified  mean  density’  p,  or  MMD  for 
short).  The  difference  between  these  two  fields  p and 
p is  again  of  second  order  in  small  quantities  and  is 
thus  very  small  where  the  TRM  theory  is  formally  valid. 
However,  the  time  derivatives  of  EMD  and  MMD  differ 
by  0(1)  amounts  because  of  the  above  discussion.  The 
MMD  is  advected  by  the  (Eulerian)  mean  flow  and  by 
the  quasi- Stokes  velocity: 

pt  + V • [(u  4-  u+)p]  = 0. 

We  shall  see  that  near  horizontal  boundaries,  the 
small-amplitude  formulae  of  McDougall  (1998)  and  MM 
to  convert  EMD  to  MMD  break  down.  In  fact,  the 
two  fields  differ  at  first,  not  second,  order  in  the  small 
quantities.  (This  is  nothing  to  do  with  the  question  of 
neutrally  stable  and  mixed  layers,  which  me  beyond  the 
scope  of  this  paper.) 

The  earliest  parameterisation  was  a simple  pair  of 
diffusivities:  + Kypzz.  This  suffers  from  the 

well-known  Veronis  effect,  in  that  fluid  within  a density 
class  is  not  conserved.  The  now  classic  GM  parameter- 
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isation  ( Gent  and  McWilliams,  1990)  takes 

Pz 

though  the  density  is  arguably  p and  not  p.  This  both 
conserves  layer  thickness  (which,  note,  is  an  integral, 
not  a conserved  quantity)  and  “smooths”  them.  It  fits 
conveniently  into  the  isopycnic  mixing  tensor  formal- 
ism of  Solomon  (1981)  and  Redi  (1982).  Visbeck  et  al. 
(1997)  give  a variant  with  horizontally  varying  diffusion 
coefficient  based  on  baroclinic  instability.  Killworth 
(1997)  uses  an  approximate  solution  to  (single  wave) 
baroclinic  instability,  obtaining  a three-dimensionally 
varying  diffusion  coefficient,  which  vanishes  for  stable 
flow.  The  scheme  mixes  q in  one  direction,  not  down- 
gradient.  Thus 

+ \ ® f^Hp\  S . 

ui'“*As;(ir/)+7'cA2 

a _ ( cos 2 8 — sin#  cos 

— sin  8 cos  9 sin2  8 ) 

and  6 is  the  angle  the  instability  wavevector  makes  with 
the  i-axis.  The  approach  is  usually  written  in  terms  of 
the  quasi-Stokes  streamfunction.  Many  other  sugges- 
tions appear  in  the  literature,  mostly  untested  ( Mar- 
shall, this  volume-MAYBE)  shows  a different  approach, 
which  is  tested  in  simple  physical  situations). 

2.2  Inferences  from  eddy-resolving  simulations 

There  have  been  several  efforts  to  use  eddy-resolving 
computations  to  enlighten  choices  about  eddy  parame- 
terisation.  To  date,  these  have  produced  not  completely 
consistent  results. 

Lee  et  al.  (1997)  indicate  that  q is  fluxed  in  a 
three-layer  channel  model,  and  not  layer  thickness. 
Roberts  and  Marshall  (2000),  however,  find  in  a depth- 
co-ordinate  model  that  the  divergent  part  of  the  tem- 
perature flux  is  not  well  correlated  with  mean  temper- 
ature gradient;  the  equivalent  for  q is  moderately  well 
correlated,  as  are  the  eddy-induced  transport  velocity. 
Wilson  (2000),  however,  in  a three-layer  channel  model 
with  forcing  varying  downstream,  finds  varying  agree- 
ment with  both  thickness  and  q flux  (thickness  having 
better  agreement  than  q)  and  large  areas  in  the  middle 
layer  in  which  the  fluxes  are  not  correlated  with  either 
gradient,  suggesting  the  importance  of  the  rotational 
component.  This  latter  is  emphasised  by  Drijfhout  and 
Hazeleger  (2001)  in  a gyral  eddy-resolving  model;  they 
show  that  the  zonally  averaged  northward  mean  thick- 
ness gradients  are  well  correlated  with  the  zonally  aver- 
aged divergent  eddy  thickness  fluxes,  while  the  equiva- 
lent for  q is  much  less  well  correlated. 


Estimates  of  eddy-induced  diffusivities  also  vary. 
There  are  strong  suggestions  that  diffusivity  varies  ver- 
tically ( Treguier  1999;  Robbins  et  al.,  2000)  as  well  as 
theoretical  suggestions  of  lateral  variation  ( Visbeck  et 
al.,  1997;  Killworth,  1997). 

The  boundary  conditions  on  surface  and  floor  of  the 
TRM  streamfunction,  which  should  apparently  be  zero 
values  in  both  locations,  also  remain  unclear.  Both 
Treguier  (1999)  and  Gille  and  Davis  (1999)  estimate 
the  streamfunction,  which  takes  extreme  values  at  both 
top  and  bottom  of  the  channels  they  considered. 

There  is  no  strong  evidence  that  any  single  glob- 
ally tested  parameterisation  (a)  gives  similar  fluxes 
(rotational  and  divergent?)  to  observed  values  from 
eddy-permitting/resolving  models  or  (b)  uniformly  im- 
proves water  masses  and  tracers.  GM  has  by  far  the 
largest  suite  of  tests,  and  while  the  evidence  is  clear 
that  - mostly  - temperature/salinity  distributions  are 
favourably  affected  by  its  use  (e.g.  Knutti  et  al.,  2000 
for  a 2.5-dimensional  model),  the  response  of  other  trac- 
ers deteriorates  ( England  and  Rahmstorf,  1999). 

This  should  not  be  surprising:  physical  eddies  have 
many  effects,  and  the  extant  parameterisations  are 
aimed  at  a small  subset  of  those  effects. 

3.  Some  details  about  TRM 
streamfunctions 

3.1  Differences  between  EMD  and  MMD 

The  TRM  approach  is  to  work  from  density  co- 
ordinates to  locate  what  “density”  variable  is  conserved 
by  a flow  consisting  of  a mean  (u)  plus  a quasi-Stokes 
velocity  (u+).  McDougall  (1998)  shows  that  this  den- 
sity (the  modified  mean  density)  is  the  inversion  of  the 
average  height  of  a density  surface, 

z(x,y,p,t)  — >•  p(x,y,z,t) 
for  some  averaging  operator.  Then 

Pt  + V-  {(u  + u+)p}  = 0. 

MM  show  that  for  small  amplitude  variability,  and 
to  second  order  accuracy, 

p = p + p = p + 0{a2) 


Pz  Pz  \PzJ 


where  the  second  term  is  usually  small  and  will  be  ne- 
glected here. 
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These  formulae  work  well  except  near  the  surface  and 
floor,  where  there  are  first-order  differences  between  p 
and  p.  These  differences  are  produced  by  the  advection 
of  fluid  laterally:  any  fluid  of  a light  density  which  is 
ever  present  at  some  (x,  y ) has  an  entry  in  the  modi- 
fied mean  density.  The  differences  occur  over  a distance 
which  is  first-order  in  perturbation  amplitude.  Fig.  1 
shows  the  near-surface  behaviour  for  a flow  whose  den- 
sity is  given  at  some  horizontal  location  by 

z - z0  = F(p  - po)  + aG(p  - po,  t) 

<5  = 0. 

(The  solution  is  for  (5(f)  = sinf.)  The  behaviour  is 
caused  by  time-averaging  near  surface  or  floor  only  be- 
ing valid  when  the  resulting  depth  z is  within  the  fluid. 

The  length  scale  in  Fig.  1 is  proportional  to  eddy 
amplitude.  For  linear  theory,  it  appears  as  a delta- 
function  boundary  layer.  When  the  eddies  have  finite 
amplitude,  the  vertical  length  scale  over  which  the  two 
densities  differ  noticeably  is  of  order 

z'  ~ p'/pz  ~a\VHp/pz\ 

which  is  the  typical  vertical  excursion  made  when  mov- 
ing a short  horizontal  distance  (a)  along  a mean  isopy- 
cnal  which  moves  significantly  vertically  only  on  the 
gyre  scale  ( L » a).  This  depth  is  small  for  ocean, 
though  not  for  the  atmosphere.  Even  with  fairly  opti- 
mistic estimates,  it  is  hard  to  produce  a vertical  scale 
much  larger  than  20  m.  So  the  distance  over  which  the 
MMD  and  EMD  differ  significantly  is  not  resolved  in 
most  climate  models , being  concentrated  in  the  last  grid 
point  (save  in  regions  such  as  the  Antarctic  Circumpolar 
Current).  Thus  the  near-boundary  differences  between 
the  two  mean  densities  will  probably  appear  to  climate 
models  as  single  grid-point  effects,  i.e.  delta  functions. 
McIntosh  and  McDougall  (1996)  show  diagnostics  from 
FRAM  which  support  this. 

No  parameterisation  of  which  I am  aware  succeeds 
in  reproducing  the  MMD  even  given  adequate  vertical 
resolution.  Fig.  2 shows  channel  model  results  for  a 
4-year  and  along-channel  average  of  an  eddy-resolving 
calculation,  with  surface  relaxation  and  parameters  de- 
signed to  permit  the  lighter  near-surface  and  denser 
near-bottom  fluid  present  in  the  MMD  to  be  resolved  by 
the  10-m  grid  spacing  (the  relevant  depth  being  50-60 
m).  Also  shown  is  a two-dimensional  calculation  using 
GM  (other  parameterisations  will  give  similar  poor  re- 
sults). It  is  clear  that  the  “pushing  forward”  of  warm 
isopycnals  present  in  the  MMD  is  not  present.  Dif- 
ferences are  very  small  at  the  lower  boundary  because 
eddy  amplitudes  were  small  there  also.  The  presence  of 


O(a)  x O(a) 


a ' (P  - Po) 

Figure  1.  The  differences  between  Eulerian  mean  and  mod- 
ified density.  The  upper  diagram  shows  that  the  densities 
are  very  close  to  each  other  in  the  fluid  interior  (differing  by 
0(a2),  where  a is  the  small  amplitude  of  the  fluctuations). 
In  a zone  of  size  a near  surface  and  floor,  the  two  densities 
differ  by  a much  larger  amount,  O(a),  as  indicated  in  the 
exploded  lower  view  (which  is  actually  the  exact  solution 
for  sinusoidal  time  variation  and  uniform  interior  density 
gradient). 


a mixed  layer  (not  treated  here)  makes  no  difference  to 
this  argument,  since  it  merely  moves  the  region  where 
EMD  and  MMD  differ  slightly  lower  (usually  to  worse 
resolution). 

Because  the  resolution  does  not  permit  the  resolution 
of  the  layer  where  MMD  and  EMD  differ,  not  all  the 
eddy  fluxes  can  be  represented.  Fig.  3 shows  the  the 
physical  near-surface  situation  contrasted  with  what  is 
represented  in  a coarse-resolution  model.  At  the  least 
density  which  the  model  is  capable  of  representing,  the 
quasi-Stokes  streamfunction  is  non-zero  (its  value,  cor- 
rectly, representing  the  lateral  eddy  flux  within  the  un- 
resolved layer,  which  thus  appears  as  a delta- function). 

An  argument  can  be  made  that  one  could  relax  the 
condition  of  = ip2  = 0 on  surface  and  floor  since  the 
model  being  run  (a)  cannot  distinguish  the  EMD  and 
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Eulerian  mean  temperature 


Modified  mean  temperature 


Figure  3.  Misrepresentation  of  quasi-Stokes  fluxes  when 
the  shallow  layer  is  not  resolved  by  the  model 


GM90  mean  temperature 


Figure  2.  The  Eulerian  and  modified  mean  density  for 
a 4-year  and  along-channel  average  of  an  eddy-permitting 
channel  model  discussed  in  the  text.  (The  average  over  the 
previous  4 year  period  is  almost  identical.)  The  problem  was 
chosen  to  provide  a larger  vertical  range  over  which  the  EMD 
and  MMD  differ  than  would  hold  for  the  real  ocean,  so  that 
the  vertical  resolution  (10  m)  was  adequate.  Also  shown 
is  a typical  two-dimensional  parameterisation  steady-state 
result,  in  this  case  following  Gent  and  McWilliams  (1990), 
using  an  eddy  diffusion  of  2000  m2  s_1.  While  the  latter 
does  not  reproduce  the  EMD  particularly  well  (true  for  a 
wide  range  of  diffusivities),  it  does  not  reproduce  the  MMD 
at  all  where  this  differs  from  the  EMD.  This  appears  to  hold 
for  most  extant  parameterisations. 


the  MMD  and  (b)  does  not  resolve  the  missing  layer. 
This  is  investigated  below. 

3.2  Mass  and  available  potential  energy 

The  differences  between  the  two  densities  have  two 
important  effects.  The  first  is  directly  concerned  with 
the  interpretation  of  mean  density.  It  is  straightforward 
to  see  that  the  low-pass  time  filtered  net  mass  in  a water 
column,  which  is  a uniquely  defined  value,  is  the  same 
whether  EMD  or  MMD  is  used: 


pdz  = / pdz  (averaging  at  constant  depth) 


= J pZp  dp  = J pzpdp  = J pdz 

(averaging  at  constant  density). 

(As  Fig.  1 suggests,  p is  lighter  at  the  surface,  but  the 
shortfall  is  made  up  at  the  floor.) 

The  same  does  not  hold  for  potential  energy,  because 
of  the  noncommutative  averaging  operators  on  products 
of  quantities.  For  small  amplitude,  the  differences  be- 
tween EMD  and  MMD  potential  energies  are  0(a2), 
and  occur  due  to  0(a2)  differences  in  the  interior  over 
a depth  range  of  order  unity,  and  O(a)  differences  over 
O(a)  depth  ranges. 
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It  is  shown  in  Killworth  (2001)  that 

,° 

APE  = / z(p  — p)dz  < 0. 
J-H 


The  difference  between  the  two  PE  expressions  lies 
in  the  variability,  fundamentally  a part  of  the  MMD.  It 
involves  an  integral  in  density  space  of  the  mean  square 
depth  fluctuations;  the  proof  is  straightforward,  either 
from  the  McDougall  (1998)  formulae  or  by  direct  eval- 
uation, and  is  not  given  here. 

Thus  the  low-pass  filtered  potential  energy  of  a fluid 
column  (a  uniquely  defined  quantity)  is  only  correctly 
evaluated  using  EMD,  and  is  consistently  underesti- 
mated using  the  MMD;  correction  terms  can  be  de- 
rived, and  involve  knowledge  of  the  variability.  This 
implies  that  energetics  cannot  consistently  be  produced 
for  coarse  models  including  parameterisations. 


4.  Creating  a parameterisation 


The  Gent  and  McWilliams  (1990)  parameterisation 
has  the  virtue  of  simplicity  and  elegance,  although  to 
function  globally  various  modifications  have  to  be  made, 
typically  those  of  ‘tapering’  neutral  density  slopes  when 
they  become  too  large.  Other  parameterisations  are 
usually  less  elegant,  and  so  need  more  modifications. 
The  depth-co-ordinate  version  of  Killworth  (1997)  is 
briefly  discussed  here  before  comparisons  are  made  in 
channel  and  global  geometries. 

The  parameterisation  uses  Unear  theory,  following 
Robinson  and  McWilliams  (1974).  Linear  theory  is  not 
always  a good  predictor  of  the  behaviour  of  a nonlinear 
eddying  system,  as  Edmon  et  al.  1980)  demonstrate, 
but  it  makes  an  indicative  starting  point  (and  holds  in 
the  same  parameter  range  as  the  McDougall  (1998)  for- 
mulae). The  system  is  assumed  to  be  slowly  varying 
in  the  horizontal,  compared  with  the  local  deformation 
radius.  The  equations  for  small  perturbations  varying 
as  exp  ik(x  cos  6 + y sin  0 - ct)  are 

(«-c)  {(;~j  -k2p^+qyP  = 0 


qy  = ficosd  — (q in  direction  6) 

u = u cos  6 + v sin  0 (velocity  in  direction  6) 


In  the  case  of  a channel  geometry,  0 is  identically 
zero. 

All  quantities  can  be  expressed  in  terms  of  the  pres- 
sure perturbation,  for  example: 


v'p' 


kcosd 

2fpo9 


Re(ipp*z);  ip2  = ~ 


kcos6 

VfilN2 


Re(ipp^) 


diffusivity  k = 


kci 


2f2P% 


_ k sin  6qy  k cos  0 qy 

Thus  to  estimate  'I',  all  that  is  needed  is  an  approxi- 
mate solution  to  the  instability  problem.  The  deforma- 
tion radius  is  estimated  from 

C=-f°  N(z)  dz,  a = C/f 
n J-H 

and  wavenumber  from  the  Eady  (1949)  result 
k = 0.51/a. 


The  angle  6 is  estimated  from  a crude  solution  max- 
imising the  growth  rate  based  on  standard  deviations 
of  u and  v over  depth.  As  shown  by  Gill  et  al.  (1974), 
6 is  small  for  east-west  flow  as  in  the  channel  case,  but 
is  near  7r/2  if  u is  small  and  v nonzero  except  in  high 
latitudes  where  /?  is  small.  The  approximate  vertical 
structure  of  k and  p is  found  by  two  iterations  of  an  ap- 
proximate (small  wavenumber)  solution  of  the  problem, 
and  k is  scaled  by  Aacz.  Here  A is  an  0(1)  constant: 
3 is  found  to  be  optimal.  The  use  of  d,  the  imaginary 
part  of  the  phase  speed,  ensures  that  there  is  no  mixing 
where  flow  is  stable,  though  there  is  always  instability 
when  v is  non-zero,  as  Pedlosky  (1987)  shows. 

The  resulting  streamfunction  possesses  two  terms. 
The  first  is  — k times  the  slope  of  the  isopycnals  nor- 
mal to  the  angle  of  instability,  which  has  similarities  to 
the  GM  formulation.  The  second,  extra  term,  has  no 
such  easy  interpretation.  Here  the  diffusivity  n varies 
three-dimensionally. 


5.  Channel  model  comparisons 

The  discussion  in  Section  3 suggests  that  the  surface 
and  floor  conditions  on  the  TRM  in  coarse  resolution 
models  are  not  obvious  (Fig.  2 being  a good  example). 
This  section  examines  solutions  to  two-dimensional  em- 
ulations of  the  three-dimensional  channel  model  of  Kill- 
worth  (1998),  using  a variety  of  formulations  to  repre- 
sent the  eddy  terms,  specifically  to  examine  the  bound- 
ary condition  question. 

The  model  covered  a longitude  range  of  2.6°,  a lat- 
itude range  of  5.2°,  centred  on  30°N,  and  a shallow 
depth  of  300  m.  The  grid  spacings  were  0.02°  east- west, 
0.018°  north-south  and  20  m vertically,  with  viscosities 
50  m2  s”1  (horizontal)  and  5 x 10~4  m2  s_1vertically 
and  diffusivities  10  m2  s-1  horizontally  and  10-4  verti- 
cally. Starting  from  a narrow  temperature  front  with 
uniform  salinity,  relaxation  towards  the  initial  temper- 
ature values  in  bands  at  north  and  south  of  the  chan- 
nel provided  a source  of  potential  energy.  This  method 


EDDY  PARAMETERISATION 


23 


has  the  advantage  that  there  are  no  regions  of  unstable 
or  neutral  stratification,  thus  avoiding  difficulties  about 
parameterisations  in  such  regions.  Averages  were  com- 
puted over  time  and  longitude  over  7.25  years  between 
days  300  and  2950. 

Two-dimensional  (latitude-depth)  simulations  were 
then  rim  using  a variety  of  two-dimensional  parameter- 
isations on  a Cartesian  grid,  and  the  4000-day  compu- 
tations (steady  in  almost  all  cases)  compared  with  the 
averages  from  the  three-dimensional  nm.  Comparisons 
were  made  with  the  temperature  field  as  a function  of  y 
(north)  and  z,  and  with  the  baroclinic  u velocity. 1 The 
comparisons  are  not  ideal.  Like  other  published  work, 
they  are  of  Eulerian  means  only,  and  over  a period  prob- 
ably an  order  of  magnitude  too  short  for  a good  statis- 
tical comparison.  (However,  the  fields  in  Fig.  3 were 
visually  unaltered  by  averaging  over  another  period  of 
similar  length,  so  the  statistics  may  be  better  than  we 
suggest.) 

Comparisons  were  made  both  visually  and  using  a 
stringent  measure  of  explained  variance  due  to  Visbeck 
et  al.  (1997).  In  each  case,  any  free  parameters  in  the 
parameterisation  were  adjusted  to  maximise  the  agree- 
ment with  the  eddy-resolving  average.  No  parameter- 
isation reproduced  the  ‘pushing  forward’  of  isopycnals 
in  the  MMD,  so  that  direct  comparisons  with  it  are  not 
useful. 

Fig.  4 shows  the  three-dimensional  time-  and  along- 
channel-averaged  solution.  To  provide  a yardstick  for 
the  various  parameterisations,  Fig.  5 shows  the  two- 
dimensional  temperature  field  using  only  advection  by 
the  actual  (two-dimensional)  velocity  fields  plus  a hori- 
zontal diffusivity  of  200  m2  s-1,  which  clearly  gives  re- 
sults very  close  to  the  three-dimensional  results. 

The  other  parameterisations  used  were  (in  order  of 
appearance  in  Figs.  6-10)  the  following: 

GM90  ( Gent  and  McWilliams , 1990,  which  has  a con- 
stant diffusivity);  Fig.  6 

K97  (more  properly,  the  depth  co-ordinate  version  of 
Killworth,  1997,  discussed  earlier,  which  computes 
a variable  diffusivity);  Fig.  7 

GMs  ( Gent  and  McWilliams , 1990,  but  with  the  stream- 
function  non-zero  at  the  surface);  Fig.  8 

Ks  ( Killworth  1997,  adapted  as  discussed  below);  Fig.9 

VP  (computing  ( v'p')y  directly  from  small-amplitude 
formulae,  also  discussed  below);  Fig.  10 


xAs  discussed  by  Killworth  (1998),  the  two-dimensional  runs 
have  no  depth-averaged  u field,  so  that  only  the  baroclinic  u can 
be  compared.  The  barotropic  u field,  as  noted  by  Killworth,  plays 
a not  inconsiderable  role  in  the  dynamics. 


Figure  4.  Contours  of  temperature  (°C;  contour  inter- 
val 0.5°  C)  for  the  time-  and  along-channel-averaged  three- 
dimensional  eddy-resolving  calculation. 


Figure  5.  Two-dimensional  simulation  of  Fig.  4,  using  a 
horizontal  diffusivity  of  200  m2  s_1 


Figure  6.  As  Fig.  4,  but  using  the  Gent  and  McWilliams 
(1990)  parameterisation,  with  k = 160  m2  s-1 
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Figure  7.  As  Fig.  4,  but  using  the  Killworth  (1997)  pa-  Figure  10.  As  Fig.  4,  but  using  a direct  parameterisation 
rameterisation  with  A - 3.  of  (v'p')y  directly  from  linear  theory,  with  A-  3 


Figure  8.  As  Fib.  4,  but  using  the  Gent  and  McWilliams 
parameterisation  modified  so  that  the  streamfunction  does 
not  vanish  at  surface  or  floor,  with  k = 1200  m2  s-1 


Latitude  (deg) 

Figure  9.  As  Fig.  4,  but  using  the  Killworth  (1997)  pa- 
rameterisation, modified  so  that  the  streamfunction  does  not 
vanish  at  surface  or  floor,  with  A = 5 


The  GM  parameterisation  is  handled  as  in  the  MOM3 
code  ( Pacanowski  and  Griffies,  1999),  with  the  TRM 
streamfunction  dropping  to  zero  over  the  last  grid  point. 
The  diffusivity  k is  taken  as  a constant.  The  K97  pa- 
rameterisation is  as  discussed  earlier.  The  third  param- 
eterisation, GMs,  attempts  to  emulate  a nonzero  value 
of  streamfunction  at  surface  and  floor.  This  is  not  an 
easy  task  numerically,  since  many  apparently  straight- 
forward approaches  generated  numerical  instabilities. 
These  included  extrapolation  of  either  the  isopycnic 
slope  or  the  streamfunction  to  the  boundary,  and  com- 
putation of  boundary  values  using  one-sided  interpola- 
tion formulae.  A slightly  unsatisfactory  approach  which 
set  the  streamfunction  at  surface  (floor)  to  the  values 
immediately  below  (above)  was  eventually  used;  the  dis- 
advantage being  that  the  v+  field  vanished  in  the  top 
and  bottom  grid  points.  The  fourth  parameterisation, 
Ks,  does  the  same  thing  for  K97. 

The  last  parameterisation,  for  the  EMD,  directly 
evaluates  {v'p')y  (very  similar  answers  are  found  for 
{v'p')y  4-  ( w'p')z ) directly  from  small-amplitude  theory, 
again  using  Killworth’ s (1997)  scalings.  Neither  calcu- 
lation requires  boundary  conditions  since  the  normal 
velocity  must  vanish  at  boundaries.  Direct  attempts 
to  parameterise  the  flux  divergence  usually  suffer  from 
Veronis  effects  (Veronis,  1975);  however,  this  approach 
does  not,  since  the  terms  have  the  same  conservation 
properties  (for  the  flux  terms)  as  the  original  system. 

The  most  accurate  version  of  the  GM90  parameter- 
isation for  this  problem  has  a k of  160  m2  s-1,  a little 
lower  than  that  cited  in  Killworth  (1998).  The  results 
for  the  GM90  (Fig.  6)  are  very  similar  to  those  of  pure 
diffusion  (5),  although  slightly  less  accurate  in  the  u 
field.  The  similarity  is  surprising  since  the  GM90  in- 
cludes the  strong  northward  (southward)  advection  near 
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the  surface  (floor)  which  is  not  present  in  the  simple  dif- 
fusive case. 

The  most  accurate  version  of  the  K97  parameterisa- 
tion  (Fig.  7)  has  A = 3,  as  used  in  Killworth  (1998)  for 
the  same  problem.  As  Fig.  7 shows,  this  parameteri- 
sation  is  the  only  one  to  produce  the  ‘doming’  of  the 
15.5°  isotherm  near  the  northern  boundary  with  any 
accuracy. 

If  w+  is  not  required  to  vanish  at  surface  and  floor, 
then  for  this  geometry  the  parameter  values  used  hith- 
erto are  insufficient  to  reproduce  the  three-dimensional 
solution.  This  is  because  the  high  northward  advec- 
tion  near-surface  is  now  lacking.  For  the  GM90  pa- 
rameterisation  (Fig.  8),  k needed  to  be  increased  an 
order  of  magnitude  (to  1200  m2  s-1)  in  order  to  repro- 
duce an  approximation  to  the  three-dimensional  fields. 
Although  the  temperature  field  looks  reasonable,  the 
corresponding  velocity  is  poorly  reproduced,  due  to  the 
strong  surface  front  near  the  southern  boundary.  A sim- 
ilar finding  holds  for  the  K97  parameterisation  (Fig.  9). 
Thus  permitting  non-zero  w+  at  surface  and  floor  has 
not  achieved  a higher  accuracy  than  maintaining  zero 
w+,  for  this  problem  and  choice  of  parameterisations. 

The  final  parameterisation  (VP)  does  not  use  the 
(v+,w+)  formulation,  but  simply  inserts  a parame- 
terisation for  mixing  directly.  The  results  (Fig.  10) 
again  show  an  accurate  representation  of  the  three- 
dimensional  result. 

In  terms,  then,  of  reproducing  the  Eulerian  mean 
density,  most  schemes  were  successful.  The  K97  and 
VP  were  only  marginally  superior  to  the  others,  and 
schemes  which  permitted  nonzero  quasi-Stokes  stream- 
functions  at  the  surface  were  quite  inferior. 


First,  there  is  also  an  element  of  tapering.  In  the  ex- 
pressions for  ipi  and  ip2  (not  shown  here)  the  first  term, 
as  noted  earlier,  includes  an  expression  which  is  the 
product  of  the  diffusivity  and  the  isopycnal  slope  nor- 
mal to  the  orientation  of  the  fastest  growing  wavenum- 
ber. This  slope  is  tapered  in  precisely  the  same  way  as 
in  isopycnal  mixing.  Second,  because  the  K97  formula- 
tion is  essentially  quadratic  in  shear  ( k increases  with 
shear,  and  operates  on  the  shear)  rather  than  linear  as 
in  GM,  an  abrupt  start,  e.g.  from  observed  tempera- 
ture and  salinity,  could  lead  to  instabilities.  The  diffu- 
sivity - and  hence  streamfunctions  - within  a column 
are  rescaled  if  necessary  so  that  k < nmax  = 104m2  s_1. 

The  role  of  baroclinic  instability  near  the  equator 
is  unclear.  While  GM  does  not  depend  on  position, 
any  scheme  involving  q mixing  must  make  choices  near 
equator  which  reflect  (a)  that  the  role  of  baroclinic  in- 
stability is  less  near  the  equator  than  at  mid-latitude 
and  (b)  that  it  is  the  east-west  density  gradient  which 
is  predominantly  creating  flows,  and  so  should  be  pref- 
erentially decreased  by  eddies.  In  addition,  preliminary 
results  analysing  OCCAM  output  (de  Vries  and  Drijf- 
hout,  personal  communication)  suggest  that  the  TRM 
streamfunction  is  much  larger  near  the  equator  than  a 
simple  constant-diffusion  GM  would  predict. 

In  any  event,  an  engineering  ‘fix’  is  required.  At 
present,  the  Coriolis  parameter  f is  taken  as  the  maxi- 
mum of  / and  some  fmin  (N.  Hemisphere  sign  conven- 
tion), and  the  deformation  radius  is  defined  from  the 
mid-latitude  and  equatorial  values  by 


amid  — C / f ; O-eq  — \J (C/2/3); 
a ——  min (timid,  Ueg)j  k — 0.51/ amid 


6.  Global  parameterisations 

Although  various  theoretical  suggestions  for  param- 
eterisations have  been  made,  relatively  few  have  been 
tested  in  global  models.  GM  has  been  shown  to  be  ro- 
bust (provided  that  tapering  arrangements  are  included 
to  prevent  poor  behaviour  in  weakly  stratified  environ- 
ments). The  K97  parameterisation  is  here  extended  to 
work  globally,  within  the  MOM3  code  (Pacanowski  and 
Griffies,  1999).  This  represents  a TRM  streamfunction 
as  an  extra  term  within  the  3x3  isopycnic  mixing 
tensor,  as 
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where  Sx,Sy  are  the  isopycnal  slopes  in  the  x and  y 
directions,  Aj  is  the  isopycnal  diffusion,  and  Ad  the 
diapycnal  diffusion. 


which  both  avoids  large  wavenumbers  and  permits  the 
approximate  solution  to  continue  to  function  reason- 
ably. 

The  last  modification  is  that  velocity  shears  are  con- 
sistently used  in  the  calculation  instead  of  density  gra- 
dients, since  thermal  wind  fails  near  the  Equator.  Near 
the  surface  these  include  the  Ekman  contribution.  This 
is  removed  in  the  surface  layer  (the  model  has  no  mixed 
layer)  by  extrapolating  the  effective  surface  velocity 
from  the  next  two  depths. 

The  above  changes  have  been  implemented  in  MOM3, 
and  run  in  a 2°  x 2°  near-global  configuration  (77°S  - 
77°N)  for  one  year  only,  since  various  features  of  the 
parameterisation  remain  under  experiment.  Monthly 
windstress  and  surface  forcing  were  employed  using 
standard  MOM3  options.  Typical  parameter  values 
were  used,  including  an  isopycnal  diffusivity  of  103 
m2  s"1.  A parallel  run  used  the  GM  formulation,  which 
is  an  option  in  MOM3,  with  k = 103  m2  s-1 . This  failed 
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after  3 months  with  erroneously  high  velocities  in  the 
Arctic.  Results  from  the  last  snapshot  files  are  shown 
from  both  runs.  Normal  diagnostics  such  as  overturn- 
ing streamfunction  are  not  useful  with  such  short  runs, 
and  for  reasons  of  space  only  one  elementary  diagnostic 
is  given  here. 

Figures  11  and  12  show  the  near-surface  (level  1)  hor- 
izontal TRM  velocities  superimposed  on  temperature 
contours  for  K97  and  GM  at  the  second  level  (depth 
37.5  m)  of  the  model,  at  the  respective  ends  of  their 
runs.  Interestingly,  neither  set  of  TRM  velocities  are 
particularly  small  compared  with  the  mean  flow.  The 
K97  TRM  flow  is  predominantly  east-west  at  latitudes 
less  than  about  40°.  This  is  caused  by  the  north-south 
orientation  of  the  fastest  growth  rate,  at  least  with  the 
approximation  used  here.  Another  approximation,  cast- 
ing the  continuous  flow  onto  a 2-layer  system  and  max- 
imising growth  rate  over  wavenumber  and  angle,  fol- 
lowed by  using  this  in  the  vertical  iteration  scheme  used 
normally  in  K97,  gives  similar  answers.  It  is  possible 
that  both  approximations  miss  other,  stronger  growth 
rates  at  higher  wavenumber.  Gill  et  al.  (1974)  give 
some  examples,  but  the  indication  from  their  work  and 
from  Pedlosky’s  (1987)  book  is  that  the  orientation  is 
expected  to  be  close  to  north-south  under  most  circum- 
stances, as  Fig.  11  suggests.  Thus  the  restriction  of  the 
TRM  flow  to  represent  q mixing  along  a single  orients 
tion  only  yields  TRM  flows  which  are  in  the  east-west 
direction  in  midlatitudes  (despite  the  fact  that  density 
gradients  are  for  the  most  part  much  stronger  north- 
south  than  east- west  in  ocean  data).  Only  in  the  subpo- 
lar areas  is  the  /3-effect  reduced  sufficiently  that  north- 
south  TRM  flows  appear,  e.g.  in  the  Antarctic  Circum- 
polar Current.  However,  the  equatorial  TRM  flows  in 
Fig.  11  are  smooth  and  well-behaved,  in  contrast  to 
the  GM  flows  in  Fig.  12  which  are  rather  noisy.  In  the 
subtropics,  both  K97  and  GM  have  weak  TRM  flows. 
However,  in  the  subpolar  region  the  GM  TRM  flows  are 
more  uniform  and  slightly  stronger  poleward  than  their 
K97  counterparts. 

7.  Discussion 

The  eddy  parameterisation  issue  is  far  from  resolved 
for  many  reasons.  Even  if  the  belief  structure  is  adopted 
that  mixing  is  done  by  eddies,  it  is  unclear  what  should 
be  mixed.  The  role  of  divergent  and  rotational  fluxes 
remains  unclear.  Coarse  models  do  not  usually  possess 
the  vertical  resolution  to  distinguish  Eulerian  and  MM 
density  near  the  surface  and  floor,  save  in  weakly  strat- 
ified regions. 

Nonetheless,  simple,  clean  parameterisations  such  as 
GM  seem  to  improve  much  (but  not  all)  of  a coarse 
model  response.  There  are  indications  that  mixing  co- 
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Figure  11.  Global  K97  results  for  near-surface  temperature 
and  horizontal  TRM  fluxes,  after  1 year’s  integration 
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Figure  12.  GM  equivalent,  after  3 months’  integration 
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efficients  vary  spatially,  but  it  is  unclear  whether  this  is 
0(1)  in  importance;  recall  that  the  channel  simulations 
showed  almost  indistinguishable  results  between  GM, 
with  a constant  diffusivity,  and  K97,  where  the  diffu- 
sivity  varied  by  two  orders  of  magnitude.  Similarly, 
more  complicated  parameterisation  schemes  appear  to 
function  as  well  as  (sometimes  better,  sometimes  not) 
simple  schemes.  Again  it  is  unclear  how  important  the 
simple-complex  distinction  is.  Certainly  it  is  not  proven 
that  any  complicated  scheme  performs  better  globally 
than  a simple  one. 

In  the  near  future,  careful  intercomparisons  of  differ- 
ent schemes  against  trusted,  statistically  reliable  three- 
dimensional  eddy-resolving  calculations  will  be  required. 
This  needs  the  community  to  invest  in  fine-resolution 
closed-gyre  forced  and  free  calculations  run  to  some 
form  of  equilibrium  to  use  as  the  testbeds  for  param- 
eterisations.  Some  of  these  need  to  be  explicitly  time- 
varying,  e.g.  on  seasonal  scales.  (Should  seasonal  varia- 
tion be  included  within,  or  without,  the  averaging  inter- 
val for  eddy  fluxes  when  applied,  say,  to  annual  mean 
coarse  resolution  modelling?)  The  community  should 
also  be  encouraged  to  seek  other  ways  of  casting  the 
parameterisation  problem  (e.g.,  the  approach  by  Mar- 
shall in  this  volume). 

Acknowledgments.  My  thanks  to  George  Nurser  for 
parts  of  Fig.  2. 

References 

de  Szoeke,  R.A.,  and  A.F.  Bennett,  Microstructure  fluxes 
across  density  surfaces.  J.  Phys.  Oceanogr.,  23,  2254- 
2264,  1993. 

Drijfhout,  S.S.,  and  W.  Hazeleger,  Eddy  mixing  of  potential 
vorticity  versus  thickness  in  an  isopycnic  ocean  model”. 
J.  Phys.  Oceanogr.  31,  481-505,  2001. 

Eady,  E.T.,  Long  waves  and  cyclone  waves.  Tellus.  1,  33-52, 
1949 

Eby,  M.,  and  G.  Holloway,  Sensitivity  of  a large-scale  ocean 
model  to  a parameterisation  of  topographic  stress.  J. 
Phys.  Oceanogr.,  24,  2577-2588,  1994. 

Edmon,  H.J.,  B.J.  Hoskins,  and  M.E.  McIntyre,  Eliassen- 
Palm  Cross  Sections  for  the  Troposphere.  J.  Atmos.  Sci., 
37,  2600-2616,  1980. 

England,  M.H.,  and  S.  Rahmstorf,  Sensitivity  of  Ventilation 
Rates  and  Radiocarbon  Uptake  to  Subgrid-Scale  Mixing 
in  Ocean  Models.  J.  Phys.  Oceanogr.,  29,  2802-2828, 
1999. 

Gent,  P.R.,  and  J.C.  McWilliams,  Isopycnal  mixing  in  ocean 
circulation  models.  J.  Phys.  Oceanogr.,  20,  150-155, 
1990. 

Gent,  P.R.,  J.  Willebrand,  T.J.  McDougall,  and  J.C. 
McWilliams,  Parameterizing  Eddy-Induced  Transports  in 
Ocean  Circulation  Models.  J.  Phys.  Oceanogr.,  25,  463- 
474,  1995. 

Gill,  A.E.,  J.S.A.  Green,  and  A.J.  Simmons,  Energy  parti- 


tion in  the  large-scale  ocean  circulation  and  the  produc- 
tion of  mid-ocean  eddies.  Deep-Sea  Res.,  21,  499-528, 
1974. 

Gille,  S.T.,  and  R.E.  Davis,  The  influence  of  mesoscale 
eddies  on  coarsely  resolved  density:  An  examination  of 
subgrid-scale  parameterisation.  J.  Phys.  Oceanogr.,  29, 
1109-1123,  1999. 

Greatbatch,  R.J.,  Exploring  the  relationship  between  eddy- 
induced  transport  velocity,  vertical  momentum  transfer, 
and  the  isopycnal  flux  of  potential  vorticity.  J.  Phys. 
Oceanogr.,  28,  422-432,  1998. 

Greatbatch,  R.J.,  and  K.G.  Lamb,  On  parameterizing  ver- 
tical mixing  of  momentum  in  non-eddy  resolving  ocean 
models.  J.  Phys.  Oceanogr.,  20,  1634-1637, 

Killworth,  P.D.,  On  the  parameterisation  of  eddy  transfer. 
Part  I:  Theory.  J.  Mar.  Res.,  55,  1171-1197,  1997. 

Killworth,  P.D.,  On  the  parameterisation  of  eddy  transfer. 
Part  II:  Tests  with  a channel  model.  J.  Mar.  Res.,  56, 
349-374,  1998. 

Killworth,  P.D.,  Boundary  conditions  on  quasi-Stokes  ve- 
locities in  parameterizations.  J.  Phys.  Oceanogr.,  31, 
1132-1155,  2001. 

Knutti,  R.,  T.F.  Stocker,  and  D.G.  Wright,  The  effects 
of  subgrid-scale  parameterizations  in  a zonally  averaged 
ocean  model.  J.  Phys.  Oceanogr.,  30,  2738-2752,  2000. 

Lee,  M.M.,  D.P.  Marshall,  and  R.G.  Williams,  On  the  eddy 
transfer  of  tracers:  advective  or  diffusive.  J.  Mar.  Res., 
55,  483-505,  1997. 

McDougall,  T.J.,  Three-dimensional  residual  mean  theory, 
pp.  269-302  in  Ocean  Modeling  and  parameterisation, 
eds.  E.  P.  Chassignet  and  J.  Verron.  Kluwer,  451  pp., 
1998. 

McDougall,  T.  J.,  and  P.C.  McIntosh,  The  temporal-residual- 
mean  velocity  Part  II:  Isopycnal  interpretation  and  the 
tracer  and  momentum  equations.  J.  Phys.  Oceanogr., 
31,  1222-1246,  2001. 

McIntosh,  P.C.,  and  T.J.  McDougall,  Isopycnal  averaging 
and  the  residual  mean  circulation.  J.  Phys.  Oceanogr., 
26,  1655-1660,  1996. 

Merryfield,  W.J.,  and  G.  Holloway,  Topographic  stress  pa- 
rameterisation in  a quasi-geostrophic  barotropic  model. 
J.  Fluid  Mech.,  341,  1-18,  1997. 

Pacanowski,  R.C.,  and  S.M.  Griffies,  MOM  3.0  manual, 
Geophysical  Fluid  Dynamics  Laboratory,  Princeton,  NJ, 
680  pp.,  1999. 

Pedlosky,  J.,  Geophysical  Fluid  Dynamics.  2nd  ed.  Springer 
Verlag,  710  pp.,  1987. 

Redi,  M.,  Oceanic  isopycnal  mixing  by  coordinate  rotation. 
J.  Phys.  Oceanogr.,  12,  1154-1157,  1982. 

Roberts,  M.J.,  and  D.P.  Marshall,  On  the  validity  of  down- 
gradient  eddy  closures  in  ocean  models.  J.  Geophys. 
Res.,  105,  28,613-28,628,  2000. 

Robbins,  P.E.,  J.F.  Price,  W.B.  Owens,  and  W.J.  Jenkins, 
The  importance  of  lateral  diffusion  for  the  ventilation  of 
the  lower  thermocline  in  the  sub-tropical  North  Atlantic. 
J.  Phys.  Oceanogr.,  30,  67-89,  2000. 

Robinson,  A.R.,  and  J.C.  McWilliams,  The  baroclinic  insta- 
bility of  the  open  ocean.  J.  Phys.  Oceanogr.,  4,  281-294, 
1974. 


28 


KILLWORTH 


Solomon,  H.,  On  the  representation  of  isentropic  mixing  in 
ocean  circulation.  J.  Phys.  Oceanogr.,  1,  233-234,  1971. 

Treguier,  A.M.,  Evaluating  eddy  mixing  coefficients  from 
eddy-resolving  ocean  models:  A case  study.  J.  Mar.  Res., 
57,  89-108.  1999. 

Treguier,  A.M.,  I.M.  Held,  and  V.D.  Larichev,  On  the  pa- 
rameterisation  of  quasi-geostrophic  eddies  in  primitive 
equation  ocean  models.  J.  Phys.  Oceanogr.,  27,  567-580, 
1997. 

Veronis,  G.,  The  role  of  models  in  tracer  studies.  In  Nu- 
merical Models  of  the  Ocean  Circulation.  Natl.  Acad,  of 
Sci.,  133-146,  1975. 

Visbeck,  M.,  J.  Marshall,  T.  Haine,  and  M.  Spall,  On  the 


specification  of  eddy  transfer  coefficients  in  coarse  reso- 
lution ocean  circulation  models.  J.  Phys.  Oceanogr.,  27, 
381-402,  1997. 

Wilson,  C.,  The  role  of  mesoscale  eddies  and  its  represen- 
tation in  numerical  models.  PhD  thesis,  University  of 
Liverpool,  128  pp.,  2000. 


This  preprint  was  prepared  with  AGU’s  IATgX  macros  v4, 
with  the  extension  package  ‘AGU++’  by  P.  W.  Daly,  version  1.6a 
from  1999/05/21. 


